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A simple model featuring a double well potential is used to represent a liquid that is quenched 
from an ergodic state into a history dependent glassy state. Issues surrounding the application of the 
Jarzynski Equality to glass formation are investigated. We demonstrate that the Jarzynski Equality 
gives the free energy difference between the initial state and the state we would obtain if the glass 
relaxed to true thermodynamic equilibrium. We derive new variations of the Jarzynski Equality 
which are relevant to the history dependent glassy state rather than the underlying equilibrium state. 
It is shown how to compute the free energy differences for the nonequilibrium history dependent 
glassy state such that it remains consistent with the standard expression for the entropy and with 
the second law inequality. 

I. INTRODUCTION 

In many real chemical and physical systems the observed distribution of the components is not that expected by 
statistical mechanics for a system at thermodynamic equilibrium. This can occur in the absence of mechanical forces 
driving the system away from equilibrium. Many compounds exist as more than one polymorph at standard room 
temperature and pressure, despite one polymorph having a significantly lower free energy than the others. Fluids can 
be cooled below the temperature at which thermodynamics predicts that a solid phase would be thermodynamically 
more stable^ Some components exist in different abundance than that predicted thermodynamically. These situations 
sometimes persist for timescales that are longer than human measurement allows or even for geological timescales. 
This behaviour might simply occur due to a slow transformation to the more stable state, or extreme rarity of 
necessary nucleating events might mean that the system is trapped in some nonergodic state for timescales that 
are incredibly long. Systems in these states are often considered as being in 'metastable' states, although it is 
perhaps inappropriate terminology for a polymorph like diamond (an allotrope of carbon that has higher free energy 
than graphite) which would normally be considered quite stable. Such systems are also often described as being in 
'nonequilibrium states', but they are nondissipative and no mechanical force is applied to prevent relaxation to the 
equilibrium distribution. To simplify terminology we will refer to these history dependent, non-dissipative, nonergodic, 
time independent nonequilibrium states as quasiequilibrium states. 

In the past it has often been assumed that subsets of the components will be equilibrated, and their relative 
distributions will be given by equilibrium, Bol tzman n distributions. Such ideas have been exploited in the so called 
energy landscape picture of the glass transition 1 2 ! 3 ! 4 ! 5 ! It has also been assumed that equilibrium thermodynamics can 
be applied to these systems, which requires that the phase space domains of the subsets do not change with small 
changes in the state point. 

Quasiequilibrium states can be formed in various ways - e.g. by temperature quenching, changing the potential 
energy function or changing the pressure of a system so rapidly that the system is not given sufficient time to adjust 
to the new conditions and the inter-domain weights can therefore not be expected to be Boltzmann. If the barrier for 
transformation between two or more local minima is high, the non-equilibrium distribution between the phase space 
domains will persist. The relative distribution between the domains will depend on the way they are prepared and 
will therefore not be a Boltzmann distribution. Nevertheless, once trapped, there is ample time for the subsystems 
to become equilibrated within their restricted phase space subdomains. Williams and Evans® produced convincing 
arguments that within these individual ergodic phase space sub-domains the internal distribution of states is given 
by a Maxwell-Boltzmann distribution and, using the fluctuation theorem as a sensitive test of aspects of the domain 
statistics, they confirmed that in aged glasses that are not too close to the glass transition, the intra-domain statistics 
are Boltzmann and the domains are robus t wit h respect to small but finite changes in the external thermodynamic 
state variables (temperature, pressure etc.) 7 ! 8 * 9 !. 

A particularly interesting example of a quasiequilibrium system is a glass. When a good glass former is prepared, 
it is not able to relax to true thermodynamic equilibrium for an extraordinarily long time, often many thousands of 
years and in the case of the natural glass obsidian, some hundred million years. The system remains in a very long 
lived, history dependent, quasiequilibrium state. Nevertheless, from a macroscopic point of view, the material appears 
to be an ordinary equilibrium solid. The fundamental thermodynamics and statistical mechanics of glass is a topic of 
active research. 
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Here we consider a simple model that could be used to represent a glass. Unlike some solids (e.g. allotropes of 
carbon), glasses have a structure which resembles that of a liquid. Because of the numerous long lived structures it is 
necessary to examine the distribution of states within the glassy system. However, the model presented below might 
also be considered to represent many other systems e.g. a protein might be frozen into a particular conformation, and 
we might be interested in the free energy of this conformation compared to the overall free energy. 

In recent years the equality for determination of free energy differences introduced by Jarzynsk i 10 * 11 ! has received 
considerable attention. This remarkable equality allows the difference in free energy between two states in thermody- 
namic equilibrium to be computed from an ensemble of nonequilibrium trajectories or pathways, of finite duration, 
which transform between the two equilibrium states. It is an interesting question as to whether this equality can 
shed light on a system which is quenched into a glassy state. It would also be interesting if this equality could be 
used to compute the difference in free energy of two polymorphs of the same compound, or to find the coexistence 
point in a phase transition. We investigate how the Jarzynski Equality can be extended to treat these systems. Since 
the Jarzynski Equality relates the free energies of different canonical states at the same temperature, we consider the 
formation of a glass by changing the potential energy surface of the system while keeping temperature constant. This 
models the formation of glassy systems by, say, altering the molecular interactions by changing the pH, or increasing 
the mole fraction of free polymer in a dispersion to form a glassy colloidal system^]. 



II. THEORY 



A. Jarzynski Equality 



Here we wi ll out line a very general approach that can be applied to arbitrary ensembles and dynamics (deterministic 
or stochastic) 13 * 14 ! It can be used to obtain the Jarzynski Equality (JE) under particular conditions, but is more 
general and will be useful in the study of the quasiequilibrium states in the next section. Consider two closed N- 
particle systems: (1) and (2) with arbitrary equilibrium distribution functions. A protocol and the corresponding 
time-dependent equations of motion are defined to transform system (1) to system (2). The dynamics may be 
thermostatted as in Eq. (|3| below or it may be isoenergetic or even unthermostatted. We define a generalised 
dimensionless "work" AA T (T;0,r), that is evaluated for a trajectory of duration r, 



cx P [AA (Dl = P %^) Zil) = /ffrodTZW 

P^(dT(r))Z^ fii\T{r))dT{T)ZW 



eq 

where is the partition function for equilibrium system i and P^q (dT) = fiq(T)dT is the probability of observing 
the infinitesimal phase volume dT, centred on the phase vector T, according to the i th equilibrium distribution 
function, /J^. The phase volume dT(r) is generated from dT using the equations of motion that take the system from 
equilibrium state (1) towards state (2) (using the forward protocol). For AX T (T) to be well defined requires that 

(a) if /£>(T) + then /i 9 2) (r(r)) ± 0, and 

(b) the converse, namely that if f$(T(r)) ^ 0, then /i^(r) ^ 0. 

The second condition is also required because if the numerator of Eq. ([TJ is zero, the "work" which is the logarithm 
of the right hand side of Eq. ([TJ will not be defined. We call these conditions the ergodic consistency conditions 
for the generalised work. They are analogous to the ergodic consistency condition for the Evans-Searles Fluctuation 
TheorerrP. 

It is trivial to prove that the exponential average of —AX r satisfies the following relation: 



(ex P (-AA)> =/ dTfW(T) f^(nr))\\dr(r)/8r\\Z^ _Z^_ 



where the brackets (. . .) x denote an ensemble average over the initial (i.e. /eg ) equilibrium distribution and 
dT(r)/dT = \\dT(T)/dT\\, This relationship is very general and shows how free energy differences can be com- 
puted from path integral information taken from nonequilibrium paths. These paths do not need to be quasistatic. 
We call this equality Eq. Q the Generalised Jarzynski Equality (GJE). 

The restriction of the integral to those regions where /jg'(r) ^ means that one completely avoids divergences in 
the function being averaged. The validity of Eq. ^ only requires: 
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(2) 

(c) an integrable region in the phase space of the final equilibrium distribution for which / eg (T(t)) ^ 0. 

We call this the ergodic consistency condition for the GJE. This condition is more general than the corresponding 
ergodic consistency condition for the generalised work. The following example illustrates a case where this condition 
breaks down and where the GJE fails. Consider the adiabatic transformation of one Hamiltonian system into a 
different Hamiltonian for which H(T(t)) ^ if (T(0)) VT. If the two equilibrium states are microcanonical and they 

(2) 

have the same energy then we will have Vr, f eq (T(t)) = 0. So in this example ergodic consistency is violated for 
both the generalised work and for the GJE. 

Of course we also make the usual physical assumptions that the dynamics is such that there are no singularities in 
the equations of motion so that the trajectories in phase space are well defined and that the Jacobian in Eq. ([5} is 
non divergent. 

To obtain the Jarzynski Equality we consider the special case of transformations using thermostatted dynamics 
between canonical equilibrium states with the same temperature. In order to determine the free energy difference, 
we consider an ensemble of initial equilibrium states at time t = that is transformed to a new state over a period 
< t < t. During this period, the ensemble of states is not at equilibrium, but if the transformation is halted at t = r, 
the system will eventually relax to a new equilibrium state. The simplest case involves a change in the functional 
form of the internal energy of the system during the period < t < t from Hq(T) to Hq 2 \t) 7 after which it is fixed 

at i?Q (r). We imagine that while these changes occur to the system of interest that it may be in contact with a very 
large heat reservoir, ensuring that the two equilibrium states are at the same temperature. If we make the system 
of interest remote from this reservoir, then it cannot possibly know the details of how the reservoir operates. As an 
example, we can model the remote reservoir by a Gaussian isokinetic reservoir where the kinetic energy of the reservoir 
particles is fixed at the value, -fQ/j e rm = 3A r t / lerm fcsT/2 where T is the equilibr ium thermodynamic temperature of 
the reservoir, N t herm is the number of thermostatted reservoir particle^ 15 ! 16 * 17 ! It is assumed that the thermostat 
temperature T is identical to the temperature of the two canonical ensembles between which we wish to calculate free 
energy differences. This reservoir can be regarded as being in thermodynamic equilibrium because it is assumed to 
have many more degrees of freedom than the system of interest. 

The equations of motion for the system during the time < t < t are written as 



Pi 

rn 



dH (q,p,X(t)) 



aSiPi 



where 

r(2) 



dH (q,p,x(t)) ^ j n j-jjggg e q ua tions A is a parametric function such that H o (\(0)) 



(3) 



and tf (A(r)) 



Hq and the function A(i) defines the transformation protocol. The switch, Si, is defined such that Si = 1 for particles 
that form the thermostatting reservoir and Si = when they are part of the system of interest. The variable a is a 
Gaussian thermostat multiplier^ that fixes the kinetic energy of the reservoir particles. It is easy to see that for such 
a system HQ herm (r,t) — ~2K t hermCt(T ,t) — Q where K t herm — ^iL±Si pf /2m is the kinetic energy of the reservoir 
particles and Q is the rate at which heat is exchanged with the synthetic thermostat. 

In this case the Liouville equation states: % = —A/ = 3iV^ erm a/, where A = ^ T is the phase space compression 
factor^ 1 ^. Hence 



0T(t) 



dr 



/e ( 9 1} (r) 
/^(r(r)) 



exp 



exp 



dtA(T(t)) 



-3N t h, 



dta(T(t)) 



exp 



(3 



dt H^"' 



(r(t)) 



(4) 



where denotes feq' evolved for a period r. In general fi 1 ' is not an equilibrium distribution. 



f (i) 



(i) 



If the equilibrium distributions fel\v) and feq'(T) are canonical and at the same temperature, it is trivial to show 
using Eq. (JTJ, that AX T /(3 is the total energy change in the system minus the energy (i.e. the heat) gained by the 
system from the thermostat (usually a negative quantity), — AQ(r; 0, r) = — J Q T dt HQ herm (T(t))]. That is, using Eqs. 
& (|4| , we see that 



f(2) 



AA r (r;0,r) = (3 / dt [H^ (T , t) 



Tjtherm t 



(r, t)} = (3 dt [flg°*(r, t)] - /3AQ T (r ; 0, r) 
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= /3 f dtH$ d {T,t)=0AW T {T;O,T). (5) 
Jo 

Here Hq d is the adiabatic (unthermostatted) time-derivative of the internal energy^. The final equality is obtained by 
consideration of the First Law of Thermodynamics, and shows that in this case (thermostatted dynamics with canonical 
initial and final distributions), AX T is just the work performed on the system in the transformation multiplied by f3: 
AX t (T;0,t) = (3AW T (T; 0, r p i i3 | i8|20 j2ii Substitution of Eq. Q into Eq. Q then gives the well known Jarzynski 
Equality, 



(exp(-/3AWg> egl = |S = exp[-(3AA], (6) 
where the partition functions are related to the Helmholtz free energy by the equation 

A = -k B T\n (^J dT exp(-/3H (T))j = -k B ThxZ. (7) 

Eq. ^ provides a way of determining the difference in the Helmholtz free energy, AA = A^ — A^, between two 
canonical equilibrium states with partition functions Z^> and by measuring the work, AW T done over a period 
r, for an ensemble of nonequilibrium pathways starting in state "1" and moving towards but not actually reaching 
equilibrium state "2". 

The same result is obtained if the initial ensemble is canonical and the dynamics is either thermostatted by a 
Nose-Hoover thermostat or the dynamics are adiabatic. For other ensembles Eq. (JlJ may not refer to work (see^). 
For example the microcanonical ensemble with the same energy Hq at times t = and t — r, the generalised "work" 
AX T , is in fact the change in heat. 

In the derivation of Eq. (|6| it is assumed that the initial distribution is given by the full canonical ensemble. The 
initial distribution must be a fully relaxed ergodic equilibrium state. The identity is then a mathematical relation 
about how the free energy difference on the right hand side is related to various integrals on the left hand side. If at 
the end of the protocol for changing the Hamiltonian the system is not in true thermodynamic equilibrium, as long as 
the ergodic consistency conditions hold, subsequent relaxation to equilibrium does not matter. If ergodic consistency 
fails (because the observed phase density at time r, has no overlap with the final equilibrium distribution) then the 
Jarzynski Equality fails. Also in any practical implementation of any GJE the generalised work needs to be properly 
defined, so the ergodic consistency condition for the generalised work takes precedence over that for the GJE itself. 

A necessary condition for Eq. ^ or Eq. (|2j to yield correct results in practise is that in the ensemble averaging 
process the time reversed path of the most probable path, must be observed. If the averaging process is not sufficiently 
exhaustive for the initial points of these possibly extremely rare events to be sampled from the initial equilibrium 
distribution, numerical evaluation of Eq. or Eq. ^ will give misleading results. One can easily see that this 
is the case. Write the work (or the generalised work) as the sum of the reversible work and the purely irreversible 
work. As the reversible work is just the free energy difference it can be taken through the average of the negative 
exponential. The average of the negative ex ponen tial of the purely irreversible work must now average to unity. This 
is just the Nonequilibrium Partition Identifying!] n { s we \\ known from the Evans Searles Fluctuation TheoremEvans 
and Searles^ that the necessary condition for this to hold in sampled data is to see the anti trajectories of the most 
probable trajectories for the process considered. This observation has an immediate impact on the calculation of free 
energy differences in the thermodynamic limit. These differences must be calculated for finite systems for a series of 
system sizes and then extrapolation must be employed in order to take the thermodynamic limit. 



B. Quasiequilibrium free energies from the quasiequilibrium partition function 

Consider an ensemble of glass-forming systems at equilibrium at t = 0. We then quench the system to a quasiequi- 
librium glass state by changing Hq over a period < t < r. After the relaxation of transients, at t = r qe > r, 
we assume that the ensemble remains in the glass state for a prolonged period of time during which the average 
properties of the system seem constant on the time scale of observation. The effect of this process on the phase space 
distribution is shown schematically in Fig. [l] If the time scale over which the system relaxes to equilibrium is very 
slow (t ^> T qe ) we may accurately model the ensemble's distribution function at t — r qe by treating it as consisting of 
a set of non-overlapping phase space domains {D a \ a = 1, N[>}. These domains partition the phase space available to 
any individual sample. By definition any sample belongs to one and only one phase space domain. Within individual 
domains the samples are ergodic (by definition) and time averaged properties are equal to ensemble averages over 
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Figure 1: A schematic diagram showing how the phase space density (represented by shading) and the location of a phase space 
volumes centred at T(0) evolve with time from an equilibrium state at time = through a period, r where the Hamiltonian is 
changing with time, and then as the system relaxes to a quasiequilibrium state at long times. 



sets of samples belonging to the same phase space domain. The domains have zero overlap - otherwise they would 
not be ergodic. When viewed separately each of these domains appears to be in equilibrium with internal weights 
given by the relative Boltzmann weights. However the relative number of ensemble members in each of these domains 
is not consistent with an equilibrium Boltzmann distribution^. Instead these relative populations are influenced by 
the details of the quench and subsequent ageing process that was used to prepare the ensemble of samples - they are 
history dependent. We call this ensemble a quasiequilibrium ensemble. The phase space distribution function for this 
ensemble has been derived by Williams and Evangel 

Following reference^, for t — r qe we can write the distribution function of a single occupied domain, D a7 as 



= s a (r)ex P [-/3ff (r)] Sa (r)ex P [-/?g (r)] = 8a (r)exp[-/?ff (r)] 

Ia[ ' J Da dTeM-PH Q (r)] JdTs a (T)eM-(3Ho(T)] Z a U 

where the switch, s a (T), is equal to unity when T E D a and zero otherwise, and 

Z a = f dTs a {T)en>[-PH (T)]. (9) 

We note that / Q (r) is the phase space density at T normalised over D a only. Only in the case of an equilibrium state 
will / Q (r) = /(r) for all T. We now write the distribution function for the quasiequilibrium ensemble, 



Safi a a (r)exp[-/?g (r)] _ j:^ iaaSa (r)ex P [-pH (r)} _ ElV a / a (r)^ 

T Nd a 

where the partition function 



E?=Wd dTexpi-PHop)] Z z Z z ' 



Z z = Y^a a drexp[-/3fr (r)]. (11) 

\ JD n 



6 



The value of a a gives the contribution of the domain a to the partition function, relative to its contribution in an 
equilibrium state. If the domain a is unoccupied, a a = 0. If we consider an equilibrium distribution of states that 
is arbitrarily partitioned into domains then we see that if the same partition function is to be obtained by summing 
over the arbitrary domains as was obtained without partitioning, then 



a a = 1, V a. 

This in turn implies the quasiequilibrium normalisation condition: 



(12) 



^2a a = N D . 

a=l 

We can define a partition function weighted free energy of an ensemble of glass samples, Az, as 



(13) 



exp[-(3A z ] = Z z = Y j a a dF exp[-/3tf (r) 

! JD„ 



Y, a <*j drs Q (r)exp[-/3iJ (r) 



^ a a Z a = y^a a exp[-(3A a ] 

a=l 



a=l 



(14) 



where the local domain free energy is A a = — ksTXxi Z a . For ease of reference we will refer to Az as the quasi- 
Helmholtz free energy. We will show later that (except at equilibrium) this free energy is not the Helmholtz free 
energy. 

Replacing the equilibrium distribution function in the definition Eq. (JTJ with the quasiequilibrium distribution 
functions, Eq. (101, and using Zz for the partition functions, it is straightforward to show that for an ensemble of 
glass samples, 



exp[AX z , T9e (r)] 



/W(r)||ar/ar(r 9B )||^ J 



(i) 



fg\T{r qe ))Z^ 



(2) 



E«=i ai 1 ] s i 1) (r )ex P [-/3g (r)]nar/^r(r ge 

(r(r ?e ))exp[-^Fo(r(r, e ))] 



n d „(2) „(2) 



7=1 



= exp [(3(H (T(r qe )) - H (T)} \\dT/dT{ Tqe ) \\ 



(2) 



(2) „(2). 



exp[(3AW T jr)] 



E^4 a) 4 2) (r(r 9e ))' 



(15) 



We note that AW T (T) = AW T (T) since no work is done during the relaxation period t < t < r qe , and AW T (r(0)) 

is the work given by Eq. In Eq. (15 I a„\ are the weights for the relaxed quasiequilibrium state which, by 
choice of r qe , are the domain weights at that same finite time. Note at this same time the intra-domain weights may 
not yet be Boltzmann. Therefore, 



AX ZiT (T)=0AW r (T) + \n 



in 



E«f4 2) (r(v)) 



(16) 
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The expression Eq. (151 was obtained from Eq. (HI and therefore AXz. Tqc (T(0)) also satisfies a modified version of 
Eq. (§: 



(exp(-AX z , T5 J) 1 = * = exp[-(3AA z ] 



,(2) 
Z_ 

7 (1) 



(17) 



Due to its resemblance to the GJE, Eq. we refer to this as the GJE for quasiequiHbrium ensembles, or simply 
the quasiequilibrium GJE. In the derivation of this relationship we have assumed: 

(d) The occupancy of the domains in the final state at t — r qe is the same occupancy as in the relaxed quasiequi- 
librium state. 

(e) For every V where f q e\T) ^ we require that fqe\T(r qe )) ^ and vice versa. Translating this into domain 



weights, if s« (r) then we require that s^\T(r qe )) ^ 0, and vice versa. 



(2) (2) 

The first assumption arises because although AW r (r(0)) does not change beyond t = r, the weights a a and s a (T(t)) 
do continue to change for t > r. This is very different to the usual circumstance for ergodic systems. We also note 
that a necessary condition for the intra-domain weights to be Boltzmann is that they relax on a time scale which is 
much shorter than that for the inter-domain weights. This implies that although the inter-domain weights do not 
significantly change for t > r qe , at the time r qe the intra-domain weights may not yet be Boltzmann. 

For an aged glass where the observable rate of change for any macroscopic property relative to the relaxation time 
of the intra-domain weights approaches infinity, we have given a proof that if the phase space domains are robust with 
respect to small changes in macroscopic parameters, the phase space distribution with in any domain is distributed in 
a Maxwell-Boltzmann distribution. Such distributions always satisfy assumption (d) above^. 

We note that if there is only one occupied domain, (e.g. a,\ = 1 and a 2 , . . . , a^ D = 0) then Az = Ai — fcsTln(l) = Ai 
where A\ is defined as in Eq. (14 1. 



If the initial and final states are both at equilibrium, AAz = AA since ai = 1V1 However, in general this is not the 
case. In order to determine the relationship between the Helmholtz free energy difference between quasiequilibrium 
states and AAz, we consider its thermodynamic definition. The Gibbs expression for the entropy is 



s = -k B J dr/(r)m[/(r)], 



(18) 



where the integral is over all space. Once the nonequilibrium process has finished, and the system slowly relaxes 
towards equilibrium (t > r) , we have every expectation that Eq. ( 18 1 obeys the second law inequality, that is 



dt ~ dt 



d(Q) 



(19) 



When the transformation process finishes no more work is done on the system and the only way the average energy 
can change is through the transfer of heat, in our case through thermostats. Given the fundamental thermodynamic 
relation for the Helmholtz free energy, 



A = (Ho) - TS, 



(20) 



and the equation for the entropy, Eq. (18 1, the Helmholtz free energy is uniquely defined. In reference^, it was shown 



that quasiequilibrium states can be treated using standard macroscopic thermodynamics. So by use of the distribution 
function, Eq. (10 1, we obtain the Helmholtz free energy of a quasiequilibrium state 



A qe = (H Q ) oe + k B T(ln{f qe (T)}) c 



(H ) qe + k B T( In 



a a s a (T)exp[-[3Ho(T)]/Z z 



(H ) qe -k B T( In 



^ s a (T)a a 

.a=l 



k B T(-/3H (T)-ln[Z z ]) 



qe 



qe 
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'Nr. 



= -k B T( In 



^ s a (T)a a 



A 2 



(21) 



where the notation (B(T)) = J D dT f qe (T)B(T) where D is all the available phase space in the glass state. Therefore, 
calculation of Az which is an ensemble average of AXz,r, e calculated along nonequilibrium trajectories, and use of 
Eq. (21 1 allows the Helmholtz free energy of a quasiequilibrium state to be determined: 



AA qe = AA Z 



k R T ( In 



" D 

E 



a (1) s (1) 



(r) 



— knT ( In 



N r 



E 



,(2)J2) 
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«V'(r(T fle )) 



(22) 



In Eq. (27) of reference [6] it was shown (for Gibbs free energies rather than Helmholtz) that the free energy 
was minimised when all the domain weights were Boltzmann, that is J D dT aa cx P(~^ g °( r )) = ,jp ex p(-/^ g o(r)) ^ 
This implies that A(ai,ci2 ■ ■ -ciNn) ls minimised when a, = 1 V i, which coincides with thermodynamic equilibrium. 
We make the standard observation of macroscopic thermodynamics that when the system is not acted on externally 
Eqs. ( |l9| & (20) give dA/dt < 0. Thus we have proved the following: as the system's distribution function moves 

V a, Eq. (19) is obeyed. By towards we mean the 
[{VA qe ■ Vg)/(Vg ■ V<?)] V<? where 
In contrast Az might 



-VA„ 



towards the equilibrium state, which is at the point a a = 1 V a, Eq. 
direction the system is moving has a component in the direction given by 
g = Yl a =i a a ~ Nri and V acts on the Nd dimensional space given by the coordinate set a 
not be a minimum in equilibrium. 



C. Quasiequilibrium free energies from a weighted sum of local domain free energies. 



A nonequilibrium free energy relation can also be obtained for systems that are quenched from an equilibrium state 
to a quasiequilibrium state by considering the average of exp(— f3AW T ) over trajectories that are in domain D a of a 
quasiequilibrium system at time r. These can then be combined to obtain a difference in free energy of the initial 
equilibrium state and the final quasiequilibrium state. In this case it is convenient to work using a weighted sum of 
distributions that are normalised over sub-domains as introduced by Williams and Evans^l 



f 9 e<X) 



E ™«/«( r ) 

Q = l 



= E 



w a s a (r)exp[-/3ff (r)] 
^ J Da dTexp[-(3H (T)} ' 



(23) 



Here Nrj is the number of domains, and w a represents the r elative w eights of these domains under the constraint 



Yla=i Wa = 1> see re f erence Ifil As above, we assume conditions (a)|(c) hold. The relationship between w a and 



be obtained by considering Eqs. (10 1 & (23 1 and is given by: 



a a f Da dTexp[-(3H (r)} _ a a Z c 
Y^cuJd dTexp[-/3ff (r)] ~ ~^ 



a a can 



(24) 



In order to develop free energy relations, we consider two possibilities: in the first case we monitor the work as the 
system is quenched from an initial equilibrium state to a quasiequilibrium state by varying Hq over a period < t < r; 
in the second case we consider the reverse process where the work is monitored as a prepared quasiequilibrium system 
is relaxed towards the ergodic equilibrium state by varying H over a period < t < r using the reverse protocol. 
In the first case we can use the relationships between the work and equilibrium canonical distribution functions, Eqs. 
([2} & ([5}, to show for any domain, D a , of a quasiequilibrium state, 



(s a (T(r qe )) ex P (-/3AVy T "(r))) 



eql 



rfi7«(r)i 



(T(r qe ))f^(T(r qe ))\\dT(r qe )/dT\\Z^ 



dT{ Tqe ) 



s a {T{T qe ))fii\T{T qe ))zW Z< 2) 

zw zw 



= exp[-/3AA 2aA ] 



(25) 
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where we have used the fact that AW^ qe (T) = AW^ qe (T) and AA 2a ,i = A^-A^\ Note that only trajectories 
that are in D a at time r qe will have a non-zero contribution to the ensemble average on the left. Using this and 
the Schrodinger-Heisenberg equivalence for phase space averages (see Section 3.3 of reference [15} we may write 
(s a (r(r ge ))) e l = (s a (T)) qe — w a . By only averaging over trajectories that are in D a at time t, we form a conditional 

ensemble average (B(r)) eq ,^ I^WM-Flr..)) Mr M)B{ r ))e 



JdTf^ q >(r)s a (r(r qc )) 



, so using Eq. (25 I we can write: 



<exp(-/3AW^ e (r))) e 



1 



Wo 



cxp[-/3AA 2a ,i]. 



(26) 



(s*(T(T qe ))) eql 

We emphasise that by (. . ) eql . a we imply that the condition is on the domain that the trajectory is in when it 
reaches the quasiequilibrium state, and not on the domain in which it starts. 



Note that in Eqs. (16 1 & (25 I, we have assumed that w a and a a do not vary with time (they remain equal to their 
value at t = r qe ). Of course this does not allow for the extremely slow relaxation to the final equilibrium state that 
might occur, on a much larger timescale. Therefore it is more accurate to say that they are constant on the accessible 
timescales. 



Conditions (a) & |(d)| specified above for the quasiequilibrium GJE are also required in this case. However, as 
discussed above, these restrictions are not likely to be of any significance for cases of practical interest, certainly not 
from the energy landscape point of view^l. 

From Eq. (26 1, and the fact that the Helmholtz free energy A — Yl a =il w o'-^ct + ksTwa ln(w a )] (see Eq. (25) of 
reference [6]) , we obtain, 



AA^ = -k B Tj2< 2) ^ (ex P (-/?AW^ e (r))) e?1;a 



(27) 



where AA eq ^ qe 



,(2) 



A^ 1 ' is the difference between the initial equilibrium Helmholtz free energy A^ and the 

(2) 

Helmholtz free energy of the nonergodic quasiequilibrium state A qe ' . 

Alternatively, we can develop an expression for the ensemble average over a single domain of the quasiequilibrium 
state, by considering the reverse process when the initial state (state (2)) is a relaxed quasiequilibrium state and 
the final state (state (1)) is an ergodic equilibrium state. Since each domain in the quasiequilibrium state is locally 
canonical, the same arguments in Section A can be used to show that for any trajectory starting at T in domain D a 
and subject to a change in Hq over a period < t < t, (rather than the longer period < t < r qe ) 



exp(pAW qe ^ eq (T)) 
The ensemble average over domain D a is given by 



fL 2 \T)\\dT/dT(T)\\Z^ 



(28) 



(ex P (-/?AWT^ e9 (r))) c 



f^(T(r))\\dT(r)/dT\\Z^ 



,(2) 



exp 



dT(r) * a (r)/W(T(r)) 
-0(AW - A^ - k B Tw a ) 



(29) 



where the Schrodinger-Heisenberg equivalence for phase space averages is used to give J dT(r) s Q (r)/ig'(r(r)) = 



(s a (T)) qe = w a and obtain the final equality. As above, we can then show that 



N D 



AA ,e^e q = _ keT ^2 W W In (exp(-/3Aiyr~ +e9 (r))) ct9e 



(30) 



In this section and Section B, we have described three approaches for determination of the free energy of the 
quasiequilibrium state that involve exponential averages of nonequilibrium path integrals, in the same way the Jarzyn- 
ski Equality is applied to equilibrium states. In the next section we will consider practical issues regarding their 
application. 
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III. MODEL AND SIMULATION DETAILS 



In order to examine the ability of the Jarzynski Equality and the new free energy expressions Eqs. (21 1 & (30 1 to 
probe the free energy of a quasiequilibrium state, we consider a simple model of a system whi ch ma y be quenched into 
a glass state. We employ a dynamical model originally developed by Hoover and coworker d 24 * 25 ! for other purposes. 
This simple dynamical system is ergodic and mixing and samples phase space canonically despite there only being 
a single particle in a one dimensional Cartesian space. To achieve this two Nose-Hoover thermostats are employed 
giving the following equations of motion 



(31) 



q = 


p_ 

rn 


P = 


m- 


Ci = 


\ m 


Cs = 


(f-- 

\ m 



-*P Z /T; 



' 3 • 



where q is the particle's position, p is its momentum, (3 — l/(/csT) where T is the average temperature regulated 
by the two thermostats and F(q) = —d$(q)/dq is the force acting on the particle. The variables t\ and T3 are time 
constants for the thermostat's feedback mechanism. 

Using the Liouville theorem, the distribution function of this system can be derived, and is given bjEl, 



t i„ n c /• \ r i r 3 exp(-/3ff e (g,p,Ci,C3))- rt9l 

/W.P.U.W (27r )3/2 (mfcBT) i/2 f^dq'e^-p^q')). 

Here H e (q,p, Ci,Cs) = H (q,p) + I^T^fCi +r|Cf ) = <S>{q) + \[p 2 /m+k B T(T 2 ( 2 + T 2 (; 2 )} where H is the Hamiltonian 
and internal energy of the unthermostatted oscillator. The partition function is 

_ (2^(mk B T)^ J^dq'eM-mq')) 2n „ 
Z e — — z, (33) 

T"lT 3 TlT 3 

where Z e is the partition function in the extended phase space. We use a double well potential to form a simple model 
of a glass. This forms a very simplistic representation of the complicated energy landscape of a real glass. It features 
a local minimum that can be separated from the global minimum on quenching. The potential, shown in Fig. [2] is 
given by the equation 

$(<?) = e (b + hq/a + b 2 q 2 /* 2 + b^/^) , (34) 

where b = 12.04541125, &i = 1.5, b 2 = —5.25 and b 4 = 0.75. We use reduced units throughout this section, where the 
length unit is c, the mass unit is m and the energy unit is k B T, resulting in the time unit a^Jm/ksT. The thermostat 
time constants have fixed values of t\ = T3 = 0.5. The potential has a global minimum of $ = at q = 1.79483214, a 
local minimum at q = —1.9385372 and a local maximum at q — 0.14370505. In Fig. |2]the potential is plotted for two 
different values of e (in units of ksT). When e = 0.1, the two local minima are separated by a barrier of the order of 

1 (i.e. ksT), a barrier that the system readily traverses. At the lower temperature, when e = 1.0, (i.e. fc^T) a system 
in the local minimum is separated from the global minimum by a barrier of the order of 6 (i.e. QksT). This is a 
significant energy barrier and the crossing of it constitutes a rare event. After quenching a large ensemble of systems 
to this low temperature, the higher energy local minimum is populated by a larger proportion of the ensemble than 
it would be when fully equilibrated. Due to the high energy barrier this situation persists for a long time. 

In the quenching experiment considered above, the initial state can be divided into two non-intersecting domains 
with q < 0.14370505 and q > 0.14370505, and both are occupied according to their equilibrium distributions: 04 = 

02 = 1. For this state Zz = Z is the usual equilibrium partition function and Az — A is the usual free energy. 
In the final state we expect that, due to our selection of the initial potential, after the quench the second domain 
(q > 0.14370505) will have an occupancy that is much higher than its Boltzmann occupancy and the first will have a 
lower than Boltzmann occupancy. The equilibrated occupancy for domain 2 is nearly zero so the change relative to 
the Boltzmann level is very large. 

Here we have deliberately selected an initial distribution that leads to 02 ^ <zi, so that we can clearly demonstrate 



the difference in information provided by the Jarzynski Equality and the new free energy expressions Eqs. (21 1 & 
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Figure 2: Potential energy function, denned by Eq. (34 l with the blue dashed line corresponding to e = 0.1 and the red solid 
line corresponding to e = 1. 



(27 1. The treatment described above can be readily applied to systems where the values of the weights for more than 



one domain are significant. 

In order to compare the Jarzynski Equality and the new relations, ensembles of 10 5 independent simulations were 
carried out starting from an initial equilibrium ensemble with e = 0.1 in Eq. (34 1. At time t — the parameter e 
was linearly increased to a final value of e = 1 at time t — r = 200. A second set of simulations was computed with 
duration of t — 2000 rather than r = 200. These times were chosen to be sufficiently short that a quasiequilibrium 
state develops, yet sufficiently long that within the domains the distribution is sufficiently close to Boltzmann at the 
end of the trajectory, that no further relaxation is required to develop the quasiequilibrium distribution. That is, 
for our model we can take r = r qe . If a more rapid protocol was used, this might not be the case and the system 
would need additional time to relax before generating the quasiequilibrium distribution. The probability distribution 
was then separated into two domains with any configuration where the position of the particle was in the range 
— oo < q < 0.14370505 being designated as in the first domain and all other configurations 0.14370505 < q < oo being 
designated as in the second domain. Recall that the value 0.14370505 is the position of the local maximum in the 
potential Eq. pi). 



IV. RESULTS AND DISCUSSION 



Distribution Functions 



The first point to be tested is whether Eq. ( 10 1 is able to accurately represent our simulations and to what degree 
they are out of equilibrium at time r = 200 with e = 1.0. In order to examine this, Fig. [3] shows f(q) obtained from 
the equilibrium distribution function, Eq. (32 I, and the distribution found in the simulations. Eq. (JTojl was then used 



to fit the data, and it was found that the best fit was obtained with the single free parameter set to ai/as = 0.03627. 
Using the normalisation condition gives a\ = 0.07000; a,2 = 1.930. 

It can be seen that Eq. (10 1 fits the data very well despite Oi and a-i being very different, which means the system is 
a long way from equilibrium by this measure. In the vicinity of the local minima in the energy (at q = 1.79483214), the 
equilibrium distribution function can be seen to have a value which is very close to zero. The observed quasiequilibrium 
distribution function has a significant value here, approximately 25 times larger than the equilibrium value. As time 
progresses, after the nonequilibrium process has finished, this difference between the two distribution functions follows 
a very slow exponential decay which will depend on the barrier height in the potential (see Fig. [2} along with the 
choice of time constants t\ and T3 in Eq. (31 1 . We will not pursue the details of this further in this paper. 

The ergodic consistency conditions for the application of the Jarzynski Equality and the quasiequilibrium G JE for 
the two quasiequilibrium states are satisfied. Firstly there is only a very slow relaxation of the weights after the 
transformations. The second condition namely that if s& (T) 7^ then we require that sf\T(T qe )) ^ 0, and vice 
versa is also satisfied. Although it is hard to see in the figure the final equilibrium distribution does have measurable 
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Figure 3: Distribution functions: the red solid line is the equilibrium distribution function with e = 1.0. The black solid circles 
are numerical data for the distribution of the quasiequilibrium state obtained by quenching the system over a period r = 200, 
and the blue dashed line is the best fit to the numerical data, obtained by adjusting the single free parameter in Eq. [To] 

density in both domains and rather more obviously in the final quasiequilibrium state there is density in both domains. 
There is also density in both domains in the initial (ergodic) equilibrium distribution function. 



B. Standard Jarzynski Equality 



By employing the standard Jarzynski Equality, Eq. we are able to compute the difference in free energy between 
the equilibrium states for e = 0.1 and e — 1.0. Of course, for the simple model under consideration, the difference in 
free energy may be readily computed using Eqs. Q, ( |33| & (34 1 for both values of e. Thus the partition function was 
evaluated numerically, and a value of AA eq — 1.7219 was obtained, where we add the subscript 'eq' to indicate that 
both states were at equilibrium. The value obtained from applying the Jarzynski Equality, Eq. (|6}, to the ensemble 
of simulations of duration r = 200, was AA = 1.7229 and for the ensemble of duration r = 2000 it was AA — 1.7231. 
This strong agreement with the equilibrium value in both cases gives a clear demonstration of how the Jarzynski 
Equality gives the free energy difference between the two equilibrated states. This might seem surprising. Despite the 
fact that the period over which the work is measured is too short to generate the final equilibrium state, the Jarzynski 
Equality refers to the state that would be reached after infinite relaxation time. Therefore, although there may be a 
long-lasting glass state, the Jarzynski Equality does not refer to that. 

Provided ergodic consistency is satisfied, the Jarzynski Equality gives the difference between the initial equilibrium 
free energy and the free energy of the final equilibrium state. However, if there are phase space domains in the final 
state that are not sampled at all, the Jarzynski Equality gives the free energy difference between the two equilibrium 
states subject to the constraint that the weights are zero in the excluded domains. The proof follows from a simple 
gedanken experiment. Suppose that the reason why those excluded domains are not sampled at all is that the potential 
energy is actually infinite for all states in those excluded domains. Then clearly the Jarzynski Equality gives the free 
energy difference between the two equilibrium states defined using the modified potential energy function for the final 
system. In this final equilibrium system all the nonzero weights are Boltzmann. 

In a real glassy state, the relaxation of the fluid is so slow that the crystalline states are never reached and rather 
the glass remains in a history dependent state for time scales beyond human experience. Applying the Jarzynski 
Equality to glass forming systems where all crystalline phases are never sampled, gives the free energy difference 
between the initial equilibrium system and the "ideal" or equilibrium glass where the weights of all the glassy phases 
are Boltzmann - both intra and inter domain weights. 

If the time over which e varied was infinitely slow, the system would be quasistatic and the process would be 
thermodynamically reversible. The amount of work, Eq. ([5}, done by every trajectory would then be the same and 
equal to the change in Helmholtz free energy and the instant the process finishes the system would be in equilibrium. 
However in the final state, when e is high, our system is not ergodic on the time scale of our simulations. By ergodic we 
mean that a single trajectory is able to sample a sufficient representation of phase space to be accurately representative 



13 




Figure 4: The distribution functions for the work done, Eq. (jHJ, for an ensemble of simulations where the quench was carried 
out over period r = 200. The distributions are shown for the work done from when the quench states for times of a) 50, b) 
100, c) 150 and d) 200. Curve e) corresponds to distribution of the work done for an ensemble of simulations when a quench 
was carried out over a period r = 2000. It shows the work done at the end of the quench (at a time of 2000). 



of the entire phase space. Recalling the measured distribution function shown in Fig. [3] this would require that e does 
not change significantly during the time it takes a single trajectory to sample a sufficient representation of the two 
peaks seen in the figure. Here the occurrence of a trajectory crossing from one peak to the other is a rare event and 
as e is increased these events become much rarer. Thus as e increases, the minimal time scale on which the system 
may not change significantly, in order to obtain something representative of a quasistatic process, diverges. 

As the process is not thermodynamically reversible it is interesting to consider the distribution functions for the 
work done by the trajectories at different times, Fig. [4] The distributions are highly skewed towards large values 
of the work by trajectories where the particle remains trapped in the local minimum of the potential through the 
quench. By comparing the distribution of the more rapidly and slowly quenched ensembles, r = 200, and r = 2000 
respectively, at the instant the quench finishes, e = 1, we gain some insight as to how the Jarzynski Equality works 
for this process. The distribution for the slower quench is much sharper and less skewed, with only a single peak, 
due to the process being significantly closer to the quasistatic limit. Thus during the slower quench many trajectories 
sample both the local minimum and the global minimum in the potential energy. In contrast to this, the more rapid 
quench is highly skewed with a second broad peak which can be observed at high values of AW T //3. This is due to 
many of the trajectories becoming stuck in either the local or global minima for prolonged times during the quench, 
i.e. loosely speaking, a break down in ergodicity. Due to the form of the exponential average in Eq. [6] the long 
skewed wing and broad second peak, for the more rapid quench, make only a small contribution to the average. This 
is exactly compensated for by the trajectories which remain trapped in the global minimum and have comparatively 
little work done on them, but make a large contribution to the average. Thus the Jarzynski Equality gives the same 
change in free energy that would be obtained from a single trajectory that is quenched quasistatically. Clearly these 
distributions are not Gaussian yet it is readily apparent that the distribution will approach a Dirac delta function as 
the quench time is extended towards infinity. 



C. Quasiequilibrium free energy expressions 



Quasiequilibrium free energy from the quasiequilibrium partition function 



As seen in Fig. [3] the distribution function of the more rapidly quenched data is fitted very well at the time of 
t = 200, with e = 1, by Eq. ^ with a 1 /a 2 = 0.03627 (o x = 0.07000; a 2 = 1.930). Using these values of a a and 
numerical integration to evaluate the local partition functions, (Z a = J dTs a (T) exp[— /3H e (T)]), Eq. (14 1 can be 
used to find that AAz = 4.276 when e is changed from e = 0.1 to e = 1 over this period. We can then compare this 



result with the value obtained using averages over nonequilibrium paths, Eq. (171 where AXz, T is given by Eq. (16 1 



This method gives AAz = 4.281, and the values are obviously in good agreement. It is worth noting that if we waited 
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Figure 5 : T he free energy for quasiequilibrium systems with a potential Eq. (34 I with e =1 and the distribution function given 
by Eq. (10 I. A minimum can be seen at equilibrium, a\ = ei2 = 1. 



long enough for the system to equilibrate and then calculated AXz, T , at this time r (which is much longer than 
200), we would once again obtain the standard change in free energy, AAz = AA eq . The large difference between 
AA eq (1.7219) and AAz (4.276) we find here shows how far the system is out of equilibrium at time r = 200 by this 
measure. 



We can now use AAz and Eq. (22 1 to determine the difference in Helmholtz free energy of the equilibrium and 



quasiequilibrium states. With a\ = 0.07000, a 2 — 1.930 and using numerical integration we find AA — 1.960, which 
is significantly different from AA eq . This shows the importance of using the new expressions for the quasiequilibrium 
free energy if it is necessary to find the free energy of the glass state. In cases where a\Z\j a 2 Z 2 — 0, it would be 
possible to obtain a good approximation to this free energy by approximating the potential energy of the glass state 
with Eq. ([34} and e = 1 for q > 0.14370505, but with $ = 00 for q < 0.14370505. However, for the data presented 
here a\Z\Ja 2 Z2 — 8.657 and this is not appropriate. 

In order to demonstrate that the Helmholtz free energy, A, given by Eq. (21 1 is minimised when a\ — a 2 — 1, as 
discussed in Section III C, we plot A as a function of a\ = 2 — a 2 for the potential, Eq. (34 1, with e = 1 in Fig. [5] 
The numerical data clearly support the theoretical result. 



2. Quasiequilibrium free energy from a weighted sum of local domain free energies 



Above we have shown that the free energy can be computed using Eqs. (17 1 & (21 1. However the approach has a 
serious drawback. If we wish to compute the free energy of a realistic model glass former using molecular dynamics 
simulations on a many body system, there will be a huge number of ergodic domains that must be considered. The 
number of domains will simply be too large to handle by these relations since it is necessary to identify what domain 
each trajectory belongs to, and to determine a and various averages involving a for each domain. A similar problem 
will occur, in general, with use of Eqs. (26 1 & (27). 



Here we devise an algorithm based on Eqs. (29 1 & (30) that avoids this. We rewrite Eq. (30) as the following 



average. 



AA qe 



*eq 



(exp(-/3AW^ e9 (r))) c 



(35) 



The over bar means that we sample master points in the phase space from the relaxed /( 2 )(r) distribution formed 
from an ensemble of quenched systems. These master points will belon g to the various domains and will by definition 



populate those domains according to the weights appearing in Eq. ( 30 1 namely Wa • For each master point we 
calculate the average (exp(— (3AW qe ^ eq (T))) a over the domain that each master point resides in. This is done by 
generating daughter points from their master. These daughter points are guaranteed to belong to the same domain 
as their master, because they are generated either by fixing the configuration from the master point and sampling 



15 



the momenta from the appropriate Maxwell-Boltzmann distribution or by simply using the equations of motion to 
ergodically generate points in that same domain. In this Monte-Carlo like procedure for correctly averaging within 
and between domains we never need to know how many domains there are, or what their explicit weights are. The 
weights occur naturally and we know how to perform averages over each ergodic subdomain. 

In the numerical work considered above, a quench time of r = 200 is sufficiently long that the numerical distribution 
is well approximated by the quasiequilibrium distribution when the quench is complete, and no additional relaxation 
time is required. This is demonstrated by the fact that the quasiequilibrium distribution function gave an excellent 



fit to the numerical data. Therefore Eq. (35 I could be applied using states generated at the end of the quench. 
The free energy difference between the quasiequilibrium and equilibrium states was determined using N q = 5 x 10 4 
trajectories or master points generated using the r = 200 quench. From each of these, an ensemble of 5000 trajectories 
was spawned to run in reverse by sampling the appropriate Gaussian distributions for the momentum p and the 



thermostat multipliers d and £3 given by Eq. (32 I. Eq. (35) was then used to calculate a change in free energy of 



AA = —AA qe ~* eq = 1.964 ±0.006. This value compares favourably with the above-mentioned directly obtained value 
of A A — 1.960, demonstrating the validity of the technique. 

V. CONCLUSIONS 

In this paper we have considered a number of different states: 

• time dependent nonequilibrium states; 

• ergodic equilibrium states where there is only one phase space domain, say domain 1, for this domain , a\ = 1 
and within this single domain phases are Boltzmann distributed; 

• nonergodic quasiequilibrium states where the domain weights a a are time independent but essentially arbitrary 
and the intra-domain weights are Boltzmann distributed; and lastly 

• constrained equilibrium states where a a = 1,0 V a , and for the occupied domains both the intra and inter domain 
weights are Boltzmann distributed. 

By studying a simple model we have shown that subject to the ergodic consistency condition, that by performing 
nonequilibrium path integrals, the Jarzynski Equality can be used to predict free energy differences between states 
that are either in thermodynamic equilibrium or constrained thermodynamic equilibrium. We have shown that subject 
to this condition this equality can even be used in systems where after the change protocol between the two states 
has completed, the final relaxation to the new equilibrium state is exceedingly slow. Our example also confirms the 
correctness of a new statistical mechanical treatmenlPof time independent, nonergodic, nondissipative nonequilibrium 
systems - so-called quasiequilibrium systems. 

We have also shown that in systems where certain phase space domains are totally unsampled in the final observed 
distribution of states (a 7 = 0), the Jarzynski Equality gives the free energy difference between the equilibrium states 
and the final constrained equilibrium state. 



We have derived three variations of the Jarzynski Equality (Eqs. (22 1, (27l & (30 1) which calculate free energy 
differences between the initial equilibrium state and the final quasiequilibrium state. Due to the intractably large 
number of domains to be considered, the first two of these new expressions would be very difficult to use on a more 



realistic examples of a glass. The third result can be restated as Eq. (35 1, which allows these problems to be overcome 
by providing the free energy of the quasiequilibrium state directly in terms of averages, without it being necessary 
to explicitly enumerate and characterise the domains. Combined, these results provide a concise illustration of how 
thermodynamics relates to glasses, polymorphs or similar arrested systems. 

The approach we have described in the main text above considers ensembles of quenched states, which in general will 
produce a number of different ergodic domains, (e.g. glasses with different physical properties or different polymorphs 
of a material). The free energies calculated therefore generally refer to the free energy of this ensemble. From a 
practical perspective, it is often of more interest to obtain the free energy of a single ergodic state (e.g. a single glass 
sample). The free energy of this system would be equal to that of a quasiequilibrium state where only one of the 



domains is populated. We note that Eqs. ( 22 1 & ( 27 1 could not strictly be applied under these conditions, since ergodic 



consistency would be violated - points in the initial equilibrium state might lead to points in the final state that are 



(2) 

not within the required domain. However, Eqs (30) & (35 I meet the required conditions, and No — 1, w\ = 1. Eq. 



(35 I takes on the particularly simple form AA a ~* eq = — fc^Tln [(exp(— /3AT / l /r "^ e<i, (r))) ce ] where a refers to the phase 



space domain that characterises the sample. This approach allows the free energy of individual ergodic subdomains 
to be determined relative to an ergodic equilibrium state. 
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